#DESCRIPTION:
#Analysis of slope limits from results of factor of safety analysis


#1. libraries and modules

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import beta,gamma
from lmfit import minimize, Parameters,report_errors
from scipy.interpolate import UnivariateSpline



#2. functions
def betadist(betaparams,B,pcF):
    #defining beta distribution parameters
    a=float(betaparams['a'].value)
    b=float(betaparams['b'].value)
    loc=float(betaparams['loc'].value)
    scale=float(betaparams['scale'].value)
    #creating fitted data
    model_pcF=beta.cdf(B,a,b,loc=loc,scale=scale)         
    #returning residual
    return (model_pcF-pcF)
    
def betadist_leastsquare_fitting(B,pcF,ai,bi,loc,scale):
    #defining power law variables
    betaparams = Parameters()
    betaparams.add('a', value=ai, vary=True, min=0.001, max=None)
    betaparams.add('b', value=bi, vary=True, min=0.001, max=None)
    betaparams.add('loc', value=loc, vary=True, min=0, max=loc)
    betaparams.add('scale', value=scale, vary=True, min=scale, max=90)
    #minimizing residuals
    result = minimize(betadist, betaparams, args=(B,pcF))
    Bfit = np.linspace(betaparams['loc'].value,betaparams['scale']+betaparams['loc'].value,90*10+1)
    pcFfit=beta.cdf(Bfit,
                    betaparams['a'].value,
                    betaparams['b'].value,
                    loc=betaparams['loc'].value,
                    scale=betaparams['scale'].value)
    
    return Bfit,pcFfit, np.sqrt(np.mean((result.residual)**2)),betaparams['a'].value,betaparams['b'].value,betaparams['loc'].value,betaparams['scale'].value


def power(powerparams,hl,B):
    hl=np.array(hl)
    B=np.array(B)
    #defining power distribution parameters
    a=float(powerparams['a'].value)
    b=float(powerparams['b'].value)
    #creating fitted data
    model_B=a*(hl**b)
    #returning residual
    return np.array((model_B-B))
    
def power_leastsquare_fitting(hl,B,ai,bi):
    #defining power law variables
    powerparams = Parameters()
    powerparams.add('a', value=ai, vary=True, min=0.001, max=None)
    powerparams.add('b', value=bi, vary=True, min=None, max=-0.00001)
    #minimizing residuals
    result = minimize(power, powerparams, args=(hl,B))
    hlfit = np.linspace(30,110,100)
    Bfit=powerparams['a'].value*(hlfit**powerparams['b'].value)
    
    return hlfit,Bfit, np.sqrt(np.mean((result.residual)**2))





def plot_results(B,pcF,a,b,loc,scale,curcolor,lstyle):
    Bfit,pcFfit,RMSE,a,b,loc,scale=betadist_leastsquare_fitting(B,pcF,a,b, loc, scale)
    plt.plot(Bfit,pcFfit*100,ls=lstyle,color=curcolor, label=hl_list[i])
    plt.axvline(x=Bfit[np.argmax(derivative(Bfit,pcFfit,1))],ls='--',lw=1,color=curcolor)
    
    #plt.axvline(x=Bfit[np.argmax(derivative(Bfit,pcFfit,2))],ls=':',lw=2,color=curcolor)
    #plt.plot(Bfit, derivative(Bfit,d_pcFfit)*100,ls='--',color=curcolor)
    #print Lith_list[j], hl_list[i], a,b,loc,scale,round(RMSE,2),round(Bfit[np.argmax(derivative(Bfit,pcFfit,1))],1), round(Bfit[np.argmax(derivative(Bfit,pcFfit,2))],1)
    return [Lith_list[j], hl_list[i], a,b, round(Bfit[np.argmax(derivative(Bfit,pcFfit,2))],1),round(Bfit[np.argmax(derivative(Bfit,pcFfit,1))],1)]
    

def derivative(x,y,order):
    s=UnivariateSpline(x,y,s=0)
    dydx=s(x,order)
    return dydx
    
    

#3. MAIN
FSdata=pd.read_csv('R32_results_summary.csv',sep=',',header=0,index_col=None)
print FSdata






hl_list=[30,50,70,90,110]
Lith_list=['CC','FC','LS', 'UM']
color_list=['b','g','r','c','m']










        
fig,ax=plt.subplots(nrows=len(Lith_list), ncols=1,sharex=True,sharey=True, figsize=(4.5,8))
main_plot=fig.add_subplot(111, frameon=False)
main_plot.set_xlabel('\n\nslope gradient, degrees')
main_plot.set_ylabel('percent failure\n\n\n')
main_plot.set_xticklabels([])
main_plot.set_yticklabels([])
main_plot.set_xticks([])
main_plot.set_yticks([])

df2=pd.DataFrame(data=None,index=None, columns=['Li','hl','a','b','f2_max', 'f1_max','10','25','50','75','100'])
for i in range(len(hl_list)):
    
    for j in range(len(Lith_list)):
               
        cur_dat=FSdata[(FSdata.hl==hl_list[i])&(FSdata.Lith==Lith_list[j])]
        B=(cur_dat['B'].values)[::-1]
        pcF=(cur_dat['pc_failure'].values)[::-1]/100.
        toe_exit=(cur_dat['toe exit'].values)[::-1]
        

        
        curax=ax[j]
        plt.sca(curax)
        curax.set_xticks((0,15,30,45,60,75,90))
        curax.set_yticks((0,25,50,75,100))
        curcolor=color_list[i]

        #initializing Beta distribution parameters for fitting
        a,b,loc,scale=1,1,0.001,89.999
        

        #print Lith_list[j]+"n"+str(hl_list[i])

        # Setting pcF=100% for slope face exits (toe exit>10)
        pcF2=np.where(toe_exit<10,pcF,np.ones(len(pcF)))
        print hl_list[i], Lith_list[j], np.round(pcF2*100,1)
        #line,=plt.plot(B,pcF2*100,'.')
        #curcolor=plt.getp(line,'color')
        dfrow=plot_results(B,pcF2,a,b,loc,scale,curcolor,'-')

        ppf=[beta.ppf(a, dfrow[2], dfrow[3], loc=0.0, scale=90.) for a in [.10,.25,.50,.75,.999]]
        dfrow={'Li':dfrow[0],'hl':dfrow[1],'f2_max':dfrow[4], 'f1_max':dfrow[5], '10':round(ppf[0]),'25':round(ppf[1]),'50':round(ppf[2],1),'75':round(ppf[3],1),'100':round(ppf[4],1), 'a':dfrow[2],'b':dfrow[3]}
        
        df2=df2.append(dfrow, ignore_index=True)

        if i==0:
            curax.set_title(Lith_list[j],fontsize='small')
        if j==len(Lith_list)-1:
            plt.legend(fontsize='x-small',loc='lower right')
        plt.ylim(-5,105)
            



df2.sort(('Li','hl'),inplace=True)
df2.reset_index(inplace=True,drop=True)
print df2
df2.to_csv('r32_pcF_analysis_results.csv',sep=',')

fig.tight_layout()


fig,ax=plt.subplots(nrows=2,ncols=1, sharex=True, sharey=True, figsize=(3.5,6))
main_plot=fig.add_subplot(111, frameon=False)
main_plot.set_ylabel('slope gradient, degrees\n\n ')
main_plot.set_xlabel('\n\nbase length, m')
main_plot.set_xticklabels([])
main_plot.set_yticklabels([])
main_plot.set_xticks([])
main_plot.set_yticks([])
for j in range(len(Lith_list)):
   
    
    cur_df2=df2[(df2.Li==Lith_list[j])]
    hl=cur_df2.hl.astype(float).values
    dydx=cur_df2.f1_max.astype(float).values
    d2ydx2=cur_df2.f2_max.astype(float).values

    #plt.plot(hl, dydx, 'ko')
    #plt.plot(hl,d2ydx2, 'kx')

    plt.sca(ax[0])

    ai,bi=30,-1

    line,=plt.plot(hl,dydx, '.')
    curcolor=plt.getp(line, 'color')
    
    
##    hl2=np.linspace(hl.min(),hl.max(),100)
##    s1=UnivariateSpline(hl,dydx,s=1,k=2)
##    dydx=s1(hl2)
##    plt.plot(hl2,dydx, ':', color=curcolor)

    hl2,dydx,RMSE=power_leastsquare_fitting(hl,cur_df2.f1_max.astype(float).values,ai,bi)
    plt.plot(hl2, dydx, '-', color=curcolor, label=Lith_list[j])
    
    if j==len(Lith_list)-1:
        plt.title('maximum pcF gradient', fontsize='small')
        plt.legend(fontsize='x-small', loc='upper right')
        ax[0].set_yticks((15,30,45,60,75))

    plt.sca(ax[1])
    plt.plot(hl,d2ydx2, '.', color=curcolor)

##    s2=UnivariateSpline(hl,d2ydx2, s=1,k=2)
##    d2ydx2=s2(hl2)
##    plt.plot(hl2,d2ydx2, ':', color=curcolor)

    hl2,d2ydx2,RMSE=power_leastsquare_fitting(hl,cur_df2.f2_max.astype(float).values,ai,bi)
    plt.plot(hl2,d2ydx2, '-', color=curcolor)

    if j==len(Lith_list)-1:
        plt.title('maximum pcF curvature', fontsize='small')
        ax[1].set_xticks((30,50,70,90,110))
        ax[1].set_yticks((15,30,45,60,75))
        ax[1].set_ylim(15,75)
        ax[1].set_xlim(25,115)

fig.tight_layout()


        
fig,ax=plt.subplots(nrows=2,ncols=2, sharex=True, sharey=True, figsize=(5,6))
main_plot=fig.add_subplot(111, frameon=False)
main_plot.set_ylabel('slope gradient, degrees\n\n  ')
main_plot.set_xlabel('\n\nbase length, m')
main_plot.set_xticklabels([])
main_plot.set_yticklabels([])
main_plot.set_xticks([])
main_plot.set_yticks([])
j2=[(0,0),(0,1),(1,0),(1,1)]
for j in range(len(Lith_list)):
   
    
    cur_df2=df2[(df2.Li==Lith_list[j])]
    hl=cur_df2.hl.astype(float).values

    pc10=cur_df2['10'].astype(float).values
    pc25=cur_df2['25'].astype(float).values
    pc50=cur_df2['50'].astype(float).values
    pc75=cur_df2['75'].astype(float).values
    pc100=cur_df2['100'].astype(float).values
    
##    hl2=np.linspace(hl.min(),hl.max(),100)
##    s10=UnivariateSpline(hl,pc10,s=1,k=2)
##    s25=UnivariateSpline(hl,pc25,s=1,k=2)
##    s50=UnivariateSpline(hl,pc50,s=1,k=2)
##    s75=UnivariateSpline(hl,pc75,s=1,k=2)
##    s100=UnivariateSpline(hl,pc100,s=1,k=2)

    hl2,s10,RMSE=power_leastsquare_fitting(hl,pc10,ai,bi)
    hl2,s25,RMSE=power_leastsquare_fitting(hl,pc25,ai,bi)
    hl2,s50,RMSE=power_leastsquare_fitting(hl,pc50,ai,bi)
    hl2,s75,RMSE=power_leastsquare_fitting(hl,pc75,ai,bi)
    hl2,s100,RMSE=power_leastsquare_fitting(hl,pc100,ai,bi)

   
    curax=ax[0,0]
    if j==len(Lith_list)-1:
        curax.set_yticks((15,30,45,60,75,90))
        curax.set_xticks((30,50,70,90,110))
        curax.set_ylim(15,90)
        curax.set_xlim(25,115)
    line,=curax.plot(hl2, s25, '-', label=Lith_list[j])
    curcolor=plt.getp(line,'color')
    curax.plot(hl, pc25, '.', color=curcolor)
    curax.set_title('25% failure',fontsize='small')
    curax.legend(loc='upper right',fontsize='x-small')


    curax=ax[0,1]
    if j==len(Lith_list)-1:
        curax.set_yticks((15,30,45,60,75,90))
        curax.set_xticks((30,50,70,90,110))
        curax.set_ylim(15,90)
        curax.set_xlim(25,115)
    line,=curax.plot(hl2, s50, '-', label=Lith_list[j])
    curcolor=plt.getp(line,'color')
    curax.plot(hl, pc50, '.', color=curcolor)
    curax.set_title('50% failure',fontsize='small')

    curax=ax[1,0]
    if j==len(Lith_list)-1:
        curax.set_yticks((15,30,45,60,75,90))
        curax.set_xticks((30,50,70,90,110))
        curax.set_ylim(15,90)
        curax.set_xlim(25,115)
    line,=curax.plot(hl2, s75, '-', label=Lith_list[j])
    curcolor=plt.getp(line,'color')
    curax.plot(hl, pc75, '.', color=curcolor)
    curax.set_title('75% failure',fontsize='small')


    curax=ax[1,1]
    if j==len(Lith_list)-1:
        curax.set_yticks((15,30,45,60,75,90))
        curax.set_xticks((30,50,70,90,110))
        curax.set_ylim(15,90)
        curax.set_xlim(25,115)
    line,=curax.plot(hl2, s100, '-', label=Lith_list[j])
    curcolor=plt.getp(line,'color')
    curax.plot(hl, pc100, '.', color=curcolor)
    curax.set_title('100% failure',fontsize='small')  


        
fig.tight_layout()
print df2
plt.show()







    
